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ABSTRACT 

In solar wind, magnetohydrodynamic (MHD) discontinuities are ubiquitous 
and often found to be at the origin of turbulence intermittency. They may also 
play a key role in the turbulence dissipation and heating of the solar wind. The 
tangential (TD) and rotational (RD) discontinuities are the two most important 
types of discontinuities. Recently, the connection between turbulence intermit¬ 
tency and proton thermodynamics has been being investigated observationally. 
Here we present numerical results from three-dimensional MHD simulation with 
pressure anisotropy and dehne new methods to identify and to distinguish TDs 
and RDs. Three statistical results obtained about the relative occurrence rates 
and heating effects are highlighted: (1) RDs tend to take up the majority of the 
discontinuities along with time; (2) the thermal states embedding TDs tend to 
be associated with extreme plasma parameters or instabilities, while RDs do not; 
(3) TDs have a higher average T as well as perpendicular temperature T^. The 
simulation shows that TDs and RDs evolve and contribute to solar wind heating 
differently. These results will inspire our understanding of the mechanisms that 
generate discontinuities and cause plasma heating. 

Subject headings: solar wind — magnetohydrodynamics — methods: numerical 
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Introduction 


The turbulent solar wind err 

ibodies discontinuit 

des fe 

Tu & Marsch 

1995; 

Marsch 

2006; 

Bruno & Carbone 

2013; 



Colburn & Sonett 

1966 


Paschmann et al. 

20R 

5). The 


tangential discontinuity (TD) and rotational discontinuity (RD) are the two most important 
yet quite different types: in the deHoffmann-Teller frame, theoretically, TDs have no normal 
components of v and B on both sides, while RDs have normal components, must obey the 
Walen relation v = ±B/^/roP, and keep |v| and |B| continuous. Furthermore, TDs allow 
jumps of density and temperature, while these parameters have the same values on both 
sides of the RDs. As to the mechanisms of how discontinuities form, there exist two kinds of 


explanations. There are 


of magnetic flux tubes flRorovskv 


empirica 


evid 


2QQSI; 


ences indica 


ung th at discontinuities are boundaries 


Miao et al. 


201 ll ). and there are suppositions that 


they form locally through n onli near interactions and 


currents (iGreco et al. 


20091). In 


Servidio et al 


may be associated with small random 


(20111) ’s 2D magnetohydrodynamic (MHD) 


simulation, most discontinuities appear to be TDs. However, in a 3D geometry, it remains 
unknown which type of discontinuity dominates. 

Dissipation of turbulence is considered an important contributor to the heating 
of the solar wind. Many rece nt studies concent rated on the role of intermittency and 


discontinuities in this process. 


Rale et al. 


along the thresholds of plasma instabilities. 


1 20091 discovere d 


Osman et al. 


stron gly enhanced fluctuations 


(20111) reported that high PVI 


(Partial Variance Increment) levels of various parameters c orrespond to i i itensiv e plasma 


heating and higher temperatures of electrons as well as ions. 


Osman et al. 


fl2012l l researched 


a large sample of data from measurements made by the Wind spacecraft and plotted the 
data distributions in the (/3||, A) parameter plane (/3|| = p\\p/A = Tj_/T]|). Thus 
they could show that the distributions are roughly bounded by curves corresponding to the 
mirror and oblique fire-hose instabilities, that the regions near the instability thresholds 
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have higher averag ed PVI, and that events with intense PVIs have A far from nnity. 


Wang et al. 


fl2013li analysed observed discontinuities with the PVI technique and found that 


the majority of them are RDs, but TDs have more obvious proton temperature increases. 
These empirical findings inspired us to investigate also the heating effects at the TDs and 
RDs obtained in our simulation. 

Numeri c al sim ulations have been employed before to understand plasma heating. 


Greco et al. 


(120081 1 assigned a path through the computational domain and then adopted 


the notion of PVI to analyse the 3D simulation data sampled along that path. Within 


their Hall MH D model they t 


intermittency. 


i ns ide ntified small-scale discontinuities being associated with 


Parashar et al. 


(1200911 demonstrated by use of a 2.5D hybrid model that an 


ion tempera ture anisotropy can be cr eated while the protons are heated by magnetic energy 


dissipation. 


Karimabadi et al 


(120131 ) conducted a full particle simulation which showed 


that, triggered by the Kelvin-Helmholtz instability with strong velocity shear, a turbulent 
cascade generates current sheets heati ng the plasma l ocally , and which yielded anisotropic 


particle distributions in that process. 


Servidio et al. 


(120141 1 allowed for a broader range of 


/3|l and the strength of the mag netic field fluc t uation s, thus obtaining results that basically 


are in accordance with those of 


Osman et al. 


(120121 1. However, neither were the data with 


high PVIs investigated, nor was the related heating of particles studied. All the mentioned 
work did also not distinguish between the various types of discontinuities in their simulation 
data. 


Motivated by all these aspects, we will here conduct a 3D numerical simulation with 
the aim to test new numerical methods to identify and analyse discontinuities, without 
assuming auxiliary paths along which data are sampled in the simulation domain. Based 
on this discontinuity identification, we present statistical results in order to investigate the 
proportion of TDs and RDs in all the discontinuities found, their parameter distribution in 
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the (/3||,^) plane, and the temperature increases in TDs and RDs. Such counts of TDs and 
RDs and their distributions have, to our knowledge, not been reported previously. These 
new simulation results will help us to understand better particle heating at intermittent 
structures in the solar wind, and thus to resolve the turbulence dissipation problem. 

In Section 2 we will describe the numerical tools and the methods employed to identify 
and categorize discontinuities. In Section 3 we present a specific case of a TD and and RD, 
as well as statistical properties of all discontinuities found in the computation domain. In 
Section 4 we shall summarize our study, and further discuss relevant physical issues. 


Methods 


2.1. Numerical Model of MHD Turbulence 


In order to evaluate the role of temperature anisotropy in MHD turbulence, we adopt 


the model which employs t 
anisotropic pressure tensor 


re compressible ic 


Meng et ah 


2012a 


eal MHD equations and incorporates an 




20131 1 : 


dp 

m 


+ V • (pv) = 0, 


dpv 

dt 


+ V • (pvv +p_lI + (P|| -p±)BB - (BB - B^\/2)/po) = 0, 


dt 


9B 

— + V X (-U X B) = 0, 
dt 

V • (p||u) + 2p||B • (B • V)u = 


5t ’ 


^ + V • (pu) + 2 p_lV • u/3 + 2(p|| - p_l)B • (B • V)u/3 = 0, 


( 1 ) 

( 2 ) 

(3) 

(4) 

(5) 


where p = (py + 2px)/3, B = B/ \B\. To describe instabilities correctly, it is c ommon to use 


a heuristic term of pre ssure relaxation restricting the temperature anisotropy flHesse fc Birn 


1992 


Birn et al. 


19951 1 denoted 5p^^/5t. In our numerical simulation we employ this 
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modified M HD model, which is not self-consistent as Vlasov models (e.g. IServidio et al 


2012, 


2QlJ). This simplification of our model, on the other hand, will help to understand 


the role of fluid-like closures of the pressure tensor, which would be dissipative with the 
pressure-transport term s. Equations ID and El keep the quantities p±/{nB) and 


adiabatically invariant (jChew et al 


19561 ). if the RHS of Equation 0] is zero. As Equation |3] 


is of the ideal MHD type, magnetic helicity is conserved. However, the model still lacks 
important kinetic physics of the solar wind, e.g. Landau damping and ion cyclotron 
resonances, which can be vital to turbulence dissipation. 


To solve the above system of equations, we utilize the BATSRUS codes (jPowell et al. 


19991: 


Toth et al. 


20121) . The simulation is conducted in a three-dimensional cartesian 


coordinate system and encompasses an absolute volume of (62.8 Mm)^ that is resolved 
in 256^ grid points. The grid resolution 250 km) is well above the ion skin depth 


100 kmb so Hal 


Rempel et al 


-dispersive physics is not included. We use the scheme proposed by 


(I2009fl in order to control numerical diffusion, and by applying such diffusion 


control strictly keep V • B = 0. This method guarantees proper dissipation and correct 
jump conditions at discontinuities. No explicit diffusive term is included in the numerical 
simulation code. Yet the magnetic and kinetic energy are subject to decay numerically. 
This decay can be attributed to the numerical scheme and grid resolution. We apply 
periodic boundary conditions to all the six surfaces of the simulation box. The initial 
conditions are set uniformly for p, T]| and T^, and randomly for v and B, with the average 
Vo = 0 and a finite guide field Bq || e^. To simulate the solar wind at 1 AU, the field 
Bo is chosen as 5 nT, while the Alfven speed is set at 50 km s“^, and p = p\\ = pi_ with 
\Jbp/ (3p) = 50 km s“^, which corresponds to a proton number density of 5 cm“^ and 
an isotropic temperature of 10^ K. For the turbulence part of th e fields we take F o urier 


amplitudes obeying the broadband initial conditions described by 


Matthaeus et al. 


fll996h . 


with |5B(k)|^ = |(5v(k)|^ oc (1 -|- fc/%nee) ^ range 10 ^m ^<A:<8xlO ^m ^ 
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The parameter Arj^^ee to be 3 x 10“^ m“^, and the spectral index q is set so that 

the slopes of the power spectral densities are both —5/3. We take -\/(?^ = 30 km s“^ 
(accordingly = 3 nT). The dimensionless cross-helicity a = 0.9. 


2.2. Methods of Selecting and Categorizing Discontinuities 


The analysis of discontinuities calls for reliable methods to select and categorize them 
in observational or numerical data. To trace abrupt spatial changes of the magnetic field, 
we use the total variance of increments (TVI) as an indicator and set a threshold for it. To 
calculate TVI, the total variance is first computed as 


AB 


(dBfiaa)\ 


( 6 ) 


where the partial derivative at grid point {i,j, k) about x is computed as 

B,{i + w,i,k) - B,{i-w,j,k) 

2w5x ■ ^ ’ 

Here 5x is the grid distance, B, denotes the corresponding component of magnetic field, 
and w is the width (in this work, we take w = 3, i.e. within a cuboid of 7^ grid points). 
For the y and derivatives, similar differences along the corresponding directions are used. 
Then the TVI is the normalized total variance 


TVI = |A5|/(|A5|), 


( 8 ) 


where the angle brackets denote the average over the whole comput ational dom a in. T his 


definition is a further development of the 


previous PV' 


(similar to the method adopted before bv iMarsch fc Tul (Il994 1). which was taken along a 


define d by 


Greco et al. 


(120081) 


given path and hence directionally sensitive. Yet the above TVI includes all directions and 
thus is unbiased for all points possibly belonging to a discontinuity. The TVI utilises more 

















information from fully 3D data, and accordingly gives more physical insight. However, 
most solar wind measurements are ID samples, so the method is inapplicable to such 
measurements. In the present work, those grid points with TVI > 3 are actually identifiable 
as discontinuity points and chosen for subsequent analysis. 


At each discontinuity point, we conduct the minimal variance analysis (MVA) 


(ISonnerup fc Cahill 


1967h in its neighbourhood, defined as a cuboid of the same given size 


for all points considered. Since Bni = is requir ed, the direction o 


can be regarded as the normal of the discontinuity llSonnerup fc Cahill 


minimal variation 


19671 ). In this work, 


we just consider a discontinuity locally as a small plane that contains the discontinuity 
point and cuts its neighbouring cuboid into two segments, so that averages of quantities on 
either side of the plane can be calculated in the segments obtained. 


To categorize the discontinuities, we use the criteria defined by 


SmithI 1 1973 ). which 


aim at judging two features: (1) whether = 0 holds, and (2) whether \B\ remains 
continuous. Hence the parameters Pi — \Bn\ /B^ and P 2 = S \B\ /B^ are calculated {B^ is 
the larger of (|S|) on both sides). The points with Pi < 0.2 and P 2 > 0.2 are categorized as 
TD, while the ones with Pi > 0.4 and P 2 < 0.2 as RD. All the other cases are categorized 
as ED(either TD or RD type, with Pi < 0.4 and P 2 < 0.2) and ND(neither TD nor RD 
type, with Pi > 0.2 and P 2 > 0.2). 


The aforementioned analysis only involves magnetic field data. To check and 
corroborate the results thus obtained, we also analyse the plasma velocity, density, and 
temperature data. In such case studies, it is trivial to check whether the density or 
temperature is continuous, but the velocity has to undergo the Walen test for RDs, i.e. one 
has to test whether v — vjjqi = ±b = ±B / y/Jiop, which is usually done statistically in a 
scatter plot. The Walen test must be conducted in the deHoffmann-Teller frame. To find its 
velocity vjj'p, in the cuboid the sum KWfc “ w) x is to be minimized. The 
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Vo tha t makes this sum a m inimum is then accepted as velocity of the deHoffmann-Teller 


frame fISonnerup et al. 


19871). 


3. Results 


The methods described above permit us to select the TDs and RDs from all cases with 
large TVI, and so we obtain corresponding data sets of discontinuities on which we can do 
individual and statistical research. 


To understand the evolution of the decaying turbulence in our MHD simulation, we plot 
in Figure [T] the squared current density {p) averaged over the whole computational domain. 
This quantity relates to the curl of the magnetic field and describes its inhomogeneity 
and energy conversion (dissipation). Its evolution in time clearly shows the following 
phases: an initial drop, s ubsequent increase and final decay. This whole trend basically 


agrees with that found by 


Matthaeus et al. 


(1199611 in their previous simulation. To further 


illustrate that process, the 2 :-component of current density is also plotted in a space-time 
display. The evolution implies that the initial drop of {p) is due to a consumption of 
the magnetic energy in the compressive turbulent plasma motion, which leads to growing 
density inhomogeneity. This increase involves the formation of thin and stretched current 
sheets (see the numerous thin and sharp structures in (b3)). Due to the decay process, the 
inhomogeneities in (b4) fade away, yet a few current-sheet-like structures remain. The trace 
power spectral density of v, b, and power spectral density of n are also shown. A region 
with spectral index close to —5/3 seems to be identified. 


To check our methods and investigate the physical properties of the discontinuities, an 
individual TD and RD are picked and illustrated in Figure |21 From the TD data we obtain 
\Bn\/\B\L = 0.17 and 5\B\/\B\l = 0.45. The TD has a normal n = (0.893, -0.355, -0.278), 








10 


almost aligned to the x-axis, and an HT velocity vht = (4.13, —10.96, —41.99) km s“^. 

In the HT frame, Panels (a) and (b) show that B and v across the discontinuity are 
quasi-parallel to the TD plane. The light grey cloud shows the sub-volume with high TVI 
and the plane is soaked therein. Panel (c) gives the TD’s Walen test, where the points are 
rather scattered, but there still exists a correlation of v and b = B/ y/nop, especially in 
the y and 2 ; components (the correlation coefficients are rxxi'i'yyiVzz = —0.06,0.88,0.82). 
For the RD we find |i?„|/|i?|L = 0.87, 5|i?|/|i?|L = 0.11, n = (0.258,-0.139,0.956) and 
vrt = (1.69,2.17,-55.56) km s“^. Panels (d) and (e) have almost identical and slightly 
bent stream lines, thus illustrating the confirmation of the Walen relation. Panel (f) 
shows the correlation coefficients {vxx^i^yy^f'zz = 0.87,0.97,0.86). The temperatures at 
the respective discontinuity points are also computed. The TD has T = 1.48 x 10^ K, 

I|l = 1.51 X 10^ K, and = 1.45 x 10® K, while the RD has T = 1.30 x 10® K, 

7|l = 1.71 X 10® K, and = 1.09 x 10® K. The TD is by 13.5% hotter than the RD in T, 
and by 33.1% hotter in T±. 

To emphasize furthermore the differences between TDs and RDs, statistical results are 
presented in Figure [3l where we plot the numbers of discontinuity points of each type as 
a function of time (left panel), as well as their percentages (right panel). At t = 0 there 
are only a few discontinuity points (53 TDs, 25 RDs, 65 EDs, 65 NDs; RDs are even not 
primary), but then the total number of discontinuity points increases with time, with RD 
becoming the dominant type. As the temporal evolution progresses, the number of TDs 
and their percentage first increase yet then decrease again slowly, behaving nearly in phase 
with that of the changing (j^), except during its initial drop phase. 

Moreover, in Figure H] we plot for TDs and RDs their beta and anisotropy locations in 
the (/3||, A) plane, where darker bins correspond to a higher number of discontinuity points, 
and with a uniform colour scale to facilitate comparisons for different times. For reference, 
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the threshold curves of the fire-hose and mirror instabilities IjHelhnger et al 


20061 1 are also 


plotted as orange and red curves, respectively. Apparently, the total numbers of TDs are 
less than those of RDs. Apart from that, they also distribute differently. At t = 60 s, the 
TD points tend to aggregate both in the centre region and near to the instability curves. 
In the decaying phase, the distribution shrinks toward the centre, and disappears finally. 
The RD points do not show this trend, with their majority still gathering around (3\\ = 0.5, 
A = 1.0, i.e. the initial values. At 60 s, some points lie beyond the instability lines but do 
not congregate there. For reference, the distributions of all grid points are supplied. They 
resemble those of the RDs. 


To investigate the heating effect of a discontinuity, the distributions of the temperatures 
found at the TD and RD points are also provided. Note that simulation of non-adiabatic 
heating is beyond the scope of this work as no real dissipation is involved. In Figure Owe 
plot the distributions of T, T\\ and at t = 60 s, with the values for TDs in red and RDs in 
blue. The TDs are slightly hotter in T and T±. The TDs have T — 1.53 x 10^ K (standard 
deviation a — 1.84 x 10^ K), whereas the RDs have T = 1.30 x 10^ K (cr = 1.68 x 10^ K). 
Also, TDs have T\_ = 1.42 x 10® K (cr = 2.37 x 10^ K), and RDs have T±_ — 1.19 x 10® K 
((T = 2.31 X 10^ K). Since both types of discontinuity evolved in the same plasma with a 
uniform initial temperature, it is reasonable to conclude that TDs tend to become hotter 
than RDs, and to infer that TDs may intrinsically be more heated than RDs. The TDs’ 
increased temperatures may be caused by plasma squeezing and adiabatic compression. In 
the present case study, the TD has a squeezing (quantified as V„n„) whose magnitude is 
5.41 times that of the RD. 





4. Conclusion 


To conclude, in this letter we describe and apply methods to identify and analyse the 
ubiquitous discontinuities appearing in 3D simulation data, and can therefore present the 
following statistical results: (1) among the identified discontinuity points, RDs represent 
the majority, (2) TDs aggregate near the extremes of the fire-hose and mirror instability 
thresholds, while RDs do not, (3) TDs are hotter than RDs, both in their T and T±. 

However, this work is also limited in some aspects. Though our simulation of decaying 
turbulence in a closed box can reveal its different stages of evolution in time, the results 
obtained are difficult to compare with turbulence observations made in the solar wind, 
where the convected plasma volume is open and can always be filled with fresh waves 
continuously supplied from a source region. To alleviate this problem (i.e. the difference 
between the reality and our case study and its temperature distributions), we selected the 
results obtained at the time when is maximal, which may represent a state being close 
to developed turbulence. 

Moreover, MHD though with thermal anisotropy lacks realistic descriptions of 
microscopic (dissipative) solar wind processes. However, the used MHD model allows 
us to reduce significantly the computational cost and to investigate discontinuities in 
three-dimensional space, a virtue of our approach which appears crucial for the analysis. 
In the identification of discontinuities, we take w = 3 and the TVI threshold as 3, values 
which are reasonably chosen but seem somehow arbitrary. Therefore, we also checked the 
outcome with a different w (ranging from 2 to 5) and TVI threshold (set lower at 2.64) 
for the statistical study, but found that the results did not change essentially, in terms of 
identification, normal direction, proportion, and distribution of all discontinuities. 

The physical mechanisms generating TDs and RDs should be further investigated. 
From another aspect, hybrid or full particle simulations should be considered, as they 
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enable the description of kinetic processes in the solar wind as well. Besides, even in 3D 
MHD simulations, the methods to study large TVI events could be further improved, e.g., 
possibly by identifying the discontinuities with model structures having a 3D geometric 
conhguration instead of by simply counting their associated points. Such approach may 
improve the identification of the discontinuities and give better statistics concerning their 
occurrence rates as well as the percentages of their local temperature increase. 


This work is supported by NSFC grants under contracts 41231069, 41174148, 
41222032, 41274132, 41474147, 41031066 and 41304133, and was carried out using the 
SWMF/BATSRUS tools developed at the University of Michigan Center for Space Environ¬ 
ment Modeling (CSEM) and made available through the NASA Comi nunitv Coordinat ed 


Modeling Center (CCMC). Figure 2 was partly produced by VAPORl Clvne et alll2007h . 


JSH, CYT, and XW are also involved in the ISSI/ISSI-BJ international team. 
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Fig. 1.— Evolution stages of the decaying turbulence in the MHD simulation, (a): Temporal 
profile of {p) in arbitrary unit, (b): Distribution of the 2 :-component of the current density 
jz in the plane z = 31.4 Mm at the given times; jz is given in arbitrary unit in agreement 
with (a), (c): Trace power spectral density of v (red) and b (blue), as well as PSD of n 
(green) at the same times. The horizontal axis represents + ky, while the vertical 

axis represents v or b(black, left) and n(green, right). Spectral indices —5/3(black) and 
—4(cyan) are also plotted in dashed lines for reference. 
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Fig. 2.— A TD (top panels) and RD (bottom panels) at t = 60 s. (a) and (b): Schematic of 
the TD, where the transparent blue plane shows the plane of discontinuity, the red and yellow 
arrows denote respectively B and v, the light grey clouds plot the TVI with thicker cloud 
representing larger values, and the red, blue and green arrows at the box corner indicate 
the X, y and 2 : directions, respectively. The background guide field is along the ^-direction, 
(c): Scatter plot of the components of v' = v — vjj'p and b = of the TD. (d) and 

(e): Schematic of the RD displayed with stream lines; cyan arrows show the direction of the 
velocity in the HT-frame. (f): Scatter plot for the RD. 
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Fig. 3.— Counts (left) of discontinuity points and their proportions (right) for each type. 
TDs, NDs, RDs and EDs are represented respectively in red, green, cyan and purple. 
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Fig. 4.— Temporal evolution of the distributions of TD points (top), RD points (centre), and 
all grid points in the whole computational region (bottom) in the (/3||,A = T^/T\\) plane, 
arranged in columns for different times. For reference, the red and orange lines give the 
thresholds for the mirror and fire-hose instabilities, respectively. The colour of a bin denotes 
the number of grid points therein. 
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Fig. 5.— Distributions of the temperatures associated with the TD and RD points. The 


left, middle and right sub-plot shows the distribution of T, Tli and T±, respectively. 







